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Abstract — A  new  noninvasive  technique  for  atrial  arrhyth¬ 
mia  analysis  is  presented  which,  based  on  time— frequency 
analysis  and  principal  decomposition,  produces  trends  of  the 
atrial  signal  characteristics  using  the  surface  ECG.  New  ho¬ 
mogeneity  measures  are  introduced  in  order  to  continuously 
measure  how  well  the  decomposed  functions  represent  the 
atrial  signals.  A  database  containing  signals  from  20  patients 
with  different  atrial  arrhythmias  (mainly  atrial  fibrillation) 
was  analyzed.  It  was  found  that  the  method  is  very  well- 
suited  for  characterizing  the  signals  and  it  is  expected  that 
the  resulting  functions  may  be  useful  for  separating  different 
arrhythmias. 

I.  Introduction 

It  is  desirable  to  find  non-invasive  methods  for  charac¬ 
terization  and  classification  of  atrial  arrhythmias,  including 
tachycardia,  flutter  and  fibrillation.  Information  contained 
in  the  atrial  activity  must,  in  some  suitable  way,  be  quanti¬ 
fied  to  accomplish  this  task.  So  far,  the  main  efforts  in  this 
field  have  been  directed  towards  atrial  fibrillation  analysis 
although  the  same  methods  in  many  cases  can  be  used  for 
flutter  and  tachycardia.  In  the  atrial  fibrillation  case,  the 
atrial  activity  in  the  ECG  have  conventionally  been  classi¬ 
fied  into  “coarse/medium/fine”  fibrillation,  see  e.g.  [1],  [2], 
to  provide  a  general  description  of  structure  and  amplitude. 
The  repetition  rate  (or  atrial  cycle  length)  of  the  f-waves  in 
the  ECG  has  also  been  investigated  and  serves  as  an  index 
of  the  degree  of  atrial  organization  [3],  [4],  Estimation  of 
the  average  repetition  rate  can  be  based  on  spectral  analy¬ 
sis.  Such  an  approach  gives  a  general  picture  of  the  signal 
by  providing  information  about  the  average  repetition  rate 
by  means  of  the  peak  location,  the  variation  in  the  rate 
by  the  width  of  the  peak  and  the  average  signal  energy  by 
the  peak  amplitude.  This  method  is  simple  but  provides 
valuable  clinical  information. 

Recently,  we  suggested  time-frequency  analysis  (TFA) 
for  a  more  detailed  temporal  characterization  of  variations 
in  the  repetition  rate,  [5].  The  potential  of  having  a  high 
temporal  resolution  was  illustrated  by  two  situations:  a  few 
patients  were  found  for  which  the  typical  irregular  rhythm 
of  atrial  fibrillation  was  suddenly  interrupted  by  short  in¬ 
tervals  of  a  different,  more  regular  rhythm.  Another  situa¬ 
tion  is  that  changes  in  repetition  rate  can  be  investigated 
during  different  interventions.  The  strategy  that  was  cho¬ 
sen  in  order  to  analyze  the  signals  on  a  second-to-second 
basis  was  to  use  an  iterative  cross-Wigner-Ville  distribu¬ 
tion  (XWVD).  In  some  patients,  the  XWVD  revealed  very 
large  variations  in  repetition  rate  (e.g.,  varying  from  5  to  7 
Hz)  whereas  very  small  variations  were  observed  in  others. 

However,  the  XWVD  models  the  frequency  variations  as 
a  frequency-modulated  sinusoid  which  has  a  low-pass  effect 


on  the  trends.  Further,  it  only  uses  the  energy  in  the  fun¬ 
damental  frequency  and  is  therefore  not  capable  of  tracking 
the  shape  of  the  signals  as  described  by  its  harmonics.  An¬ 
other  limitation  is  that  the  computational  complexity  is 
relatively  high. 

Atrial  signals  may  be  nonstationary  but  are  repetitive 
and  thus  they  can  during  short  intervals  be  represented  by 
a  fundamental  which  reflects  the  repetition  rate  and  a  har¬ 
monic  pattern  which  reflects  the  shape  of  the  hbrillatory 
waveform.  One  way  to  achieve  a  detailed  feature  extrac¬ 
tion  in  the  time-frequency  plane  for  this  type  of  signals  is 
to  decompose  the  time-frequency  distribution  into  descrip¬ 
tive  functions  (’’principal  components”):  spectral  profile, 
frequency  shift  trend  and  amplitude  scaling  trend,  [6]. 

This  paper  addresses  the  problem  of  characterizing  dif¬ 
ferent  atrial  arrhythmias  using  a  modified  TFA-based  prin¬ 
cipal  component  decomposition  of  residual  ECGs  in  order 
to  extract  clinically  interesting  features  such  that  different 
arrhythmias  or  different  intervals  during  one  arrhythmia 
can  be  differentiated.  The  goal  is  to  find  a  set  of  trends 
that  can  serve  as  a  basis  for  classification  of  atrial  arrhyth¬ 
mias. 

II.  Methods 

We  propose  a  TFD-basecl  characterization  of  residual 
ECGs.  The  methods  consists  of  four  parts:  a  preprocessing 
stage  which  performs  prefiltering  and  QRST  cancellation  in 
order  to  produce  an  atrial  signal.  In  the  second  stage,  the 
time-frequency  distribution  is  computed  and  decomposed 
into  a  set  of  descriptive  functions.  Finally,  a  postprocess¬ 
ing  stage  is  used  to  analyze  the  spectral  content  of  the 
frequency  trend. 

A.  Preprocessing 

A  spatiotemporal  QRST  cancellation  scheme  was  used 
[7].  Since  the  spectral  content  of  interest  in  the  resulting 
residual  ECG  signal  is  well  below  25  Hz,  the  residual  ECG 
can  be  downsampled  from  1  kHz  to  50  Hz.  This  operation 
considerably  reduces  the  amount  of  data  to  be  processed. 
The  main  parts  of  the  ventricular  activity  (QRS  complex 
and  T  wave)  is  removed  in  the  residual  ECG  signal  and  thus 
this  signal  contains  primarily  atrial  activity  i.e.  P  waves 
during  sinus  rhythm  and  f-waves  during  atrial  fibrillation. 

B.  Time-frequency  analysis 

In  the  present  paper,  we  focus  on  atrial  tachycardia, 
atrial  flutter  and  atrial  fibrillation  and  for  this  purpose  only 
the  frequency  interval  between  2.5-25  Hz  is  of  interest.  This 
interval  is  chosen  such  that  slow  tachycardias  down  to  2.5 
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in  common  that  they  are  repetitive.  During  tachycardia  the 
cycles  are  often  of  equal  length.  For  fibrillation,  the  cycle 
length  may  vary  very  rapidly  and  the  signal  is  referred  to 
as  being  nonstationary.  However,  for  short  intervals  (typ¬ 
ically  1-2  seconds)  atrial  fibrillation  signals  can  be  viewed 
as  approximately  stationary.  The  result  is  that  when  an¬ 
alyzing  the  spectral  content  of  all  these  types  of  different 
atrial  signals  on  a  second-to-second  basis  the  result  will  be 
a  fundamental  and  a  number  of  harmonics  reflecting  the 
shape  of  the  waveforms.  The  harmonic  pattern  have  negli¬ 
gible  energy  above  25  Hz.  Another  reason  to  use  the  lower 
limit  at  2.5  Hz  is  that  possible  small  QRST  residuals  ap¬ 
pear  at  low  frequencies  depending  on  the  RR  interval  and 
the  window  length.  If  the  QRST  residuals  are  of  high  am¬ 
plitude  or  occur  in  a  periodic  way  they  may  have  spectral 
content  above  2.5  Hz. 

Let  the  frequency  vector  f  =  \fo  ■■■fNf-i]T  denote  an 
increasing  sequence  (length  Nf)  of  normalized  frequencies. 
The  corresponding  DFT  matrix,  F,  (Nf-by-Nw  where  Nw 
is  the  signal  window  length)  can  then  be  written 

F  =  [  1  e~^2nt  e~^27Tt2  ...  j  ^ 

where  1  is  a  column  vector  of  length  Nf.  If  the  frequencies 
are  uniformly  distributed  the  DFT  matrix  is  unitary.  In  the 
present  method,  the  Short-Term  Fourier  transform  (STFT) 
with  a  logarithmic  frequency  scale  for  the  frequency  in¬ 
terval  2.5-25  Hz  such  that  a  doubling  in  frequency  for  all 
frequencies  corresponds  to  the  same  number  of  frequency 
bins.  The  motivation  for  this  is  to  be  able  to  compare  the 
harmonic  pattern  for  intervals  with  different  fundamental 
frequency  (or  atrial  cycle  length).  A  signal  window  length 
of  Nw  =  128,  i.e,  about  2.5  seconds,  is  used.  The  frequency 
vector  was  defined  as 

i  Nf~ 1 

f  =  [  2.5  2.5  •  1077/  ...  2.5-lO^R  ]  (2) 

where  Nf  =  Nw.  Further,  the  diagonal  entries  in  the  di¬ 
agonal  matrix  H  (Nw-by-Nw)  represent  a  window  function 
(here  chosen  as  Hamming).  The  observed  signal,  x(n),  is 
represented  by  a  data  matrix,  X,  in  which  the  columns  are 
overlapping  data  sequences.  For  a  window  distance  L  and 
a  total  signal  length  of  K  window  intervals  (total  signal 
length  is  thus  (K  —  1)  —  by  —  L  +  Nw)  the  data  matrix 
(Nw-by-K)  is  given  by 

;c(0)  .  .  .  X((K  ~  1)£) 

X  =  :  :  (3) 

.  x(Nw  -  1)  ...  x((K  -  1  )L  +  NW-1)  . 

where  column  k  represents  the  signal  in  the  k: th  window 
interval. 

The  STFT  of  the  observed  signal  x(n)  can  be  written  as 

Q  =  FHX  (4) 

where  Q  is  Nf-by-K  and  the  fc:th  column  contains  the 
spectrum,  qfc,  for  the  fc:th  window  interval. 

(5) 


A  time-frequency  distribution,  Q,  of  a  nonstationary  but 
cyclic  signal  can,  using  principal  decomposition,  be  divided 
into  three  trends:  a  spectral  profile  function,  0,  a  frequency- 
shift  function,  9k,  and  an  amplitude  function,  ak,  using  an 
iteration  method,  which  shifts  the  spectra  for  the  different 
window  intervals  such  that  the  first  principal  component 
of  the  correlation  matrix  based  on  all  spectra  represents 
as  much  of  the  energy  as  possible,  [6].  This  is  achieved 
when  the  peaks  for  all  spectra  are  at  the  same  location. 
The  shift  needed  for  each  of  the  spectra  is  the  frequency- 
shift  function,  Ok-  The  inner  product  between  the  first 
principal  component  and  each  spectrum,  shifted  such  that 
the  peaks  match,  reflects  the  amplitude,  a k,  of  the  signal  in 
that  window  interval.  The  iterative  procedure  is  initialized 
with,  Q1-0)  =  Q ,i  =  0.  The  correlation  matrix,  Rb),  for 
the  distribution  ,  Qb),  is  given  by 

Rb)  =  Qb)Q  (i)H  (6) 

The  eigenvector  corresponding  to  the  largest  eigenvalue, 
0b),  can  be  calculated  using  the  iterative  ’’power  method”. 
An  extended  eigenvector,  0b),  with  0  (maximum  fre¬ 
quency  shift)  extra  samples  in  both  the  beginning  and  the 
end  is  created  in  order  to  allow  selection  of  different  parts  of 
the  vector.  These  extra  samples  are  set  to  the  same  value 
as  the  start  and  end  samples  of  0b)  respectively.  There 
are  now  20  +  1  possible  sets  of  Nf  samples  that  can  be 
selected  from  0b).  The  selection  is  done  using  a  frequency 
shift  matrix  Jg  defined  as 

Je  =  [OiV/xie+e)  I  NfxNf  OW/X(e-0)]  (7) 

which  selects  Nf  samples  from  Nf  +  20. 

Maximization  of  the  scalar  product  between  each  spec¬ 
trum,  in  and  all  possible  shifts  of  the  eigenvector, 
Jg0b),  with  respect  to  9  is  done  using  a  grid  search  of  9 
in  the  interval  [—0,0]  in  order  to  find  the  frequency  shift 

0b)  =  argmax|(J^b))Hq(o)|  (8) 

U 

A  (i) 

The  corresponding  amplitude  scaling  function  af  is  then 
the  inner  product 

df  =  (V,0b))Hqbb  (9) 

k 

The  vector  q^  is  prolonged  with  0  extra  samples  in  both 
the  beginning  and  the  end  in  the  same  way  as  for  the  eigen¬ 
vector  above  and  each  shifted  spectrum,  qb+J\  in  Qb+J) 
is  constructed  as 

qii+i)  =  J_e-oqi0)  (io) 

k 

The  purpose  of  this  shift  is  to  generate  a  TFD  with  a  higher 
energy  concentration  around  the  dominant  frequency  in  or¬ 
der  to  calculate  an  even  more  concentrated  spectral  profile 
(principal  eigenvector  of  the  correlation  matrix). 


Q=[q0  ...  qjv-i] 
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to  detect  convergence  of  the  procedure.  When  the  differ¬ 
ence  in  total  energy  between  the  present  and  previous  iter¬ 
ation  is  less  than  10%  of  the  total  energy  for  the  previous 
iteration,  the  iterations  are  terminated.  Otherwise  i  =  i+ 1 
and  the  procedure  is  repeated  from  eqn.  (6). 

At  the  final  iteration,  a  homogeneity  function ,  <5fc ,  is  cal¬ 
culated  in  order  to  measure  how  well  the  spectral  profile 
represents  each  spectra  in  the  time-frequency  distribution 

4  ”  II  (W(l))  llll  q[0)  II  (11) 

k 

This  measure  is  calculated  both  for  the  entire  frequency 
interval  (2.5-25  Hz)  representing  mainly  the  fundamental 
frequency  (dqfc)  but  also  for  the  frequency  interval  9-25  Hz 
which  contains  the  harmonics  in  order  to  measure  how  well 
the  shape  matches  the  shape  represented  by  the  spectral 
profile  (S-2,k)- 

D.  Postprocessing 

One  example  of  postprocessing  is  to  investigate  modu¬ 
latory  properties  in  the  atrial  signals.  This  can  be  done 
by  performing  spectral  analysis  of  the  frequency  trend,  §k, 
i.e.  the  frequency  of  the  fundamental  peak  in  the  spectral 
profile,  <f>  =  [4>o  ■  ■  ■  4>Nf-i],  adjusted  by  the  frequency  shift 
function 

^  -^argmaXj(/>i-4  (12) 

The  ”  modulation”  spectrum  S  of  §k  is  calculated  using  the 
FFT. 


III.  Results 

In  order  to  investigate  the  properties  of  different  atrial 
arrhythmias,  the  proposed  method  was  applied  in  20 
recordings  from  patients  with  atrial  arrhythmias.  Three 
patients  had  atrial  tachycardia  and  17  had  chronic  atrial 
fibrillation.  Three  patients  were  studied  during  rhythm 
controlled  respiration  (0.125  Hz)  for  the  purpose  of  inves¬ 
tigating  possible  modulatory  effects  on  the  fibrillation  fre¬ 
quency.  Each  recording  was  of  one  minute  duration  and 
lead  V\  was  investigated  (although  V2  and  V3  were  also 
used  for  QRST  cancellation). 

The  performance  of  the  proposed  method  is  demon¬ 
strated  by  four  examples  which  are  representative  for  the 
entire  database.  Five  seconds  from  each  of  the  four  ECG 
signals  in  which  the  ventricular  activity  has  been  cancelled 
are  shown  in  Fig.  1.  The  corresponding  decomposed  time- 
frequency  distributions  are  presented  in  Figs.  2-5.  Each  of 
these  figures  shows  the  logarithmic  STFT  in  the  upper  left 
corner.  The  resulting  spectral  profile,  (f>,  is  plotted  below 
the  STFT;  for  comparison  the  power  spectrum  is  also  pre¬ 
sented.  In  the  upper  right  corner  the  three  trends  9k,  a,k 
and  Sk  are  plotted.  In  the  rightmost  plot  the  two  homo¬ 
geneity  functions  5itk  and  S2,k  are  plotted.  Finally,  in  the 
lower  right  corner  the  spectrum,  S,  of  the  frequency  trend, 
Sfc,  is  shown  and  is  of  particular  interest  in  examples  3 
and  4. 
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Fig.  1.  Five  seconds  from  each  of  the  four  ECG  signals  in  which  the 
ventricular  activity  has  been  cancelled. 


A.  Example  1 

Figure  2  shows  a  typical  case  of  atrial  fibrillation  with 
a  dominant  fundamental  frequency  at  around  6  Hz.  The 
frequency  shift  function  reveals  a  relatively  large  frequency 
variation.  In  this  case,  two  harmonics  are  visible  in  the 
spectral  profile.  From  the  homogeneity  functions  in  the 
rightmost  plot  it  is  noted  that  the  spectral  profile  relatively 
well  represents  each  spectra  in  the  distribution  both  con¬ 
cerning  the  fundamental  and  the  harmonic  pattern  (the 
average  values  of  8\^  and  S2y k  are  0.94  and  0.90  respec¬ 
tively).  This  means  that  the  shape  of  the  fibrillation  waves 
are  relatively  constant  independently  of  frequency  and  am¬ 
plitude.  In  this  example  no  obvious  patterns  can  be  seen  in 
the  amplitude  function  and  the  frequency  trend  spectrum. 
Comparing  the  power  spectrum  of  the  entire  one-minute 
signal  to  the  spectral  profile,  it  is  evident  that  the  funda¬ 
mental  peak  of  the  latter  spectrum  is  narrower  and  that 
the  harmonics  can  be  more  easily  discerned. 
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Fig.  2.  Example  1:  Decomposition  of  a  1  minute  atrial  fibrillation 
time- frequency  representation.  Upper  right:  frequency  shift  func¬ 
tion,  amplitude  function  and  homogeneity  function  (2.5-25  Hz  - 
solid  line,  9-25  Hz  -  dotted  line).  Lower  left:  Spectral  profile 
-  solid  line,  power  spectrum  -  dotted  line.  Lower  right:  Power 
spectrum  of  the  frequency  trend  Sfc. 


B.  Example  2 

The  second  example  is  shown  in  Fig.  3.  In  this  case  the 
typical  irregular  rhythm  of  atrial  fibrillation  was  suddenly 
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can  be  easily  detected  in  both  the  constant  value  of  the 
frequency  shift  trend  but  also  in  the  sudden  increase  of  the 
amplitude  trend.  Similar  to  example  1,  two  harmonics  are 
observed  but  the  peaks  are  here  more  distinct.  Again,  a 
large  difference  is  found  between  the  power  spectrum  and 
the  spectral  profile.  It  is  interesting  to  note  that  both  ho¬ 
mogeneity  measures  are  slightly  closer  to  one  during  the 
regular  rhythm.  In  average  the  homogeneity  values  were 
Si  =  0.92  and  Si  =  0.88  indicating  that  the  shape  given 
by  the  spectral  profile  also  is  representative  during  weaker 
irregular  rhythm.  No  obvious  peak  is  detected  in  the  fre¬ 
quency  trend  spectrum.  However,  it  is  curious  to  note  that 
the  largest  peak  above  0.1  Hz  occur  at  the  normal  respi¬ 
ratory  rate  (0.31  Hz  and  0.32  Hz  respectively)  in  Figs.  2 
and  3. 
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Fig.  3.  Example  2:  Decomposition  of  a  1  minute  atrial  fibrillation 
time-frequency  representation  with  a  10  s  long  regular  rhythm  at 
30  s.  (See  legend  to  Figure  1  for  further  explanations.) 
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Fig.  4.  Example  3:  Decomposition  of  a  1  minute  atrial  fibrillation 
time-frequency  representation  during  rhythm  control  respiration 
(0.125  Hz).  (See  legend  to  Figure  1  for  further  explanations.) 
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Fig.  5.  Example  4:  Decomposition  of  a  1  minute  atrial  tachycardia 
time-frequency  representation  during  rhythm  control  respiration 
(0.125  Hz).  (See  legend  to  Figure  1  for  further  explanations.) 


C.  Examples  3  and  4 

Figure  4  shows  a  decomposed  atrial  fibrillation  signal 
with  two  harmonics  in  the  spectral  profile  and  with  rel¬ 
atively  large  frequency  shifts  (0k)-  In  Fig.  5,  the  same 
decomposition  is  performed  for  an  atrial  tachycardia  case 
(3  harmonics  combined  with  an  almost  constant  frequency 
shift  function  (0k))-  A  number  of  observations  can  be 
made:  First,  the  homogeneity  measures  are  as  expected  in 
average  very  close  to  one  (i>i  =  0.99  and  Si  =  0.97)  for  the 
more  regular  arrhythmia  in  Fig.  5,  but  also  for  the  atrial 
fibrillation  case  in  Fig.  4  these  measures  are  relatively  high 
(ji  =  0.95  and  Si  =  0.91).  Secondly,  in  both  cases  obvi¬ 
ous  peaks  are  present  around  the  respiratory  frequency  at 
0.125  Hz  in  the  frequency  trend  spectra.  This  peaks  dis¬ 
appeared  when  blocking  the  autonomous  nervous  system 
with  atropine.  Finally,  when  the  rhythm  is  stationary  the 
power  spectrum  and  the  spectral  profile  are  identical. 

IV.  Conclusions 

Analysis  of  signals  from  a  database  of  20  patients  with 
different  arrhythmias  (mainly  atrial  fibrillation)  shows  that 
the  proposed  method  is  very  well-suited  for  characterizing 
the  atrial  signals  as  indicated  by  the  homogeneity  trends. 


Comparing  the  power  spectrum  of  the  signals  to  the  spec¬ 
tral  profiles,  it  is  evident  that  the  peaks  of  the  spectral 
profiles  are  narrower  and  that  the  peaks  of  the  harmon¬ 
ics,  representing  atrial  morphology,  can  be  much  better 
discerned.  Further,  modulatory  properties  in  the  atrial  ac¬ 
tivity  can  be  obtained  from  the  frequency  trend. 
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